Perturbative evolution of particle orbits around Kerr black holes: time domain 

calculation 
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We consider the problem of the gravitational waves produced by a particle of negligible mass 
orbiting a Kerr black hole. We treat the Teukolsky perturbation equation in the time domain nu- 
merically as a 2+1 partial differential equation. We model the particle by smearing the singularities 
in the source term by the use of narrow Gaussian distributions. We have been able to reproduce ear- 
lier results for equatorial circular orbits that were computed using the frequency domain formalism. 
The time domain approach is however geared for a more general evolution, for instance of nearly 
geodesic orbits under the effects of radiation reaction. 
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I. INTRODUCTION 



It will not be long before the current set of interferometric gravitational wave detectors such as LIGO, VIRGO, 
GEO and TAMA achieve sensitivities at which detections are expected to be possible with a reasonable event rate. 
The ground-based interferometric detectors are now starting full-blown "science runs" and the plans for more sensitive 
second and even third generation detectors are in the drawing boards. Even LISA, the space-based observatory, seems 
■ to be a possibility in a reasonably near future, 
l/^ \ This work tries to help in setting the stage for realistic calculations of what now seems to be a very likely occurrence 
in active galactic nuclei: the capture of compact stellar-sized objects by the supermassive black holes lurking in the 
center of most galaxies. This type of events will generate gravitational waves within LISA's planned sensitivity band. 

Since the captured objects are small compared to the capturing black hole, they can be treated in a first approxi- 
mation by black hole perturbation theory. But one thorny problem still remains that has stalled progress in this area. 
We can model the waves coming from many kinds of stationary particle orbits around the black hole, but real orbits 
will not be stationary. The emission of gravitational waves will back-react on the object making the orbit decay in 
Y \ time. One of the challenges in this area is to find workable prescriptions for these radiation reaction forces that can be 
used in numerical simulations that can produce accurate waveforms and energy fluxes for observers at large distances 
from the hole to be used as templates for the interferometric detectors. The prescriptions that have been worked 
J> . out so far involve calculating quite complex integrals along the past light cone of the particle 0, Q] • Since almost all 
perturbative treatments of particles around rotating holes are done in the frequency domain, applying these radiation 
reaction schemes has proved challenging. Frequency domain computations become computationally expensive when 
one is dealing with realistic non-circular orbits. 

We try here to work out a simple case of equatorial circular particle orbits around a rapidly rotating black hole to 
show how to treat this problem perturbatively in the time domain. We hope this might help to open the door to the 
use of the proposed fully relativistic radiation reaction schemes in the calculation of the more complex but interesting 
problem of modeling the orbital decay and gravitational radiation emission of particles in more general elliptic and/or 
inclined orbits. Up to now, almost all treatments have been based on the separability of the equation in the frequency 
domain and calculate the energy and waveforms for the first few £ multipoles of the spheroidal harmonics expansion 
once the radial part has been dealt with. 

Using the frequency domain approach, quite a few detailed simulations of the gravitational waves emitted by a 
particle in a bound orbit around a black hole have been carried out in recent years. The first results for radiation 
emitted by a particle in orbit around a Schwarzschild black hole were carried out by Detweiler, who pioneered these 
techniques [sj- In a series of six papers, Poisson and various collaborators studied in detail the gravitational wave 
emission to infinity and into the horizon of a particle in circular orbit around a non-rotating hole 0, 0, @, 0; S 13 • 

Other groups have been doing extensive work on Post-Newtonian expansions for various types of particle orbits 
around Kerr holes (see |lfj|| and references therein). Work has been done on eccentric orbits around Schwarszchild 
|ll| , circular orbits slightly inclined away from the equatorial plane of Kerr holes |12| , and particles with spin orbiting 
Kerr holes on equatorial circular orbits |13|. 

For single static orbits (ignoring the decay by radiation reaction) to get good accuracies in the frequency domain 
approach involves summing for values of I up to 12 for quite a few of the above-harmonic frequencies, requiring 
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Authors 


Type of result 


Detweiler (1978) 
Poisson (1993,1995) 


equatorial 
circular orbits around Kerr BH 


Kojima & Nakamura (1984) 


infall trajectories with 
ang mom 7^ 


Sasaki & Nakamura (1982) 


infall along symmetry 
axis of Kerr BH 


Tanaka et al (1996) 


spinning particles in 
circ equatorial orbits 


Hughes (2000) 


slightly eccentric and 
non-equatorial orbits 



TABLE I: A summary of some of the more relevant frequency domain results fora binary black hole system in the particle limit 
that have appeared in the recent literature. 

close to 3000 iterations using a numerical algorithm that involve complex Green's function calculations [T^j . Further 
discussion of the number of i modes needed, and also a discussion of the self-force in the case of a scalar field can be 
seen in the paper by Burko |l5j . 



II. TIME DOMAIN EVOLUTION OF THE TEUKOLSKY EQUATION 



We will perform the numerical evolution with the latest version of a code originally developed by Laguna and 
collaborators 0] to evolve the Tcukolsky equation in the time domain from a generic initial data slice. The method 
analyzes the radiation at "infinity" by dealing with the s = —2 version of the equation. Then the equation is rewritten 
as a first order matrix equation and numerically integrated. The code has been tested treating scalar fields, scattering 
of gravitational waves and analysis of quasi-normal ringing and power law tails of the outgoing radiation. 

We start with the original Teukolsky equation ^3] before separation is performed, 
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+ [s 2 cot 2 9-s]ifj = 4vr(r 2 + a 2 cos 2 9)T 
" 2 — 2Mr + a 2 postulating an ansatz for the tp function of the form, 

ip = *(r,t, 6) e"™*. 

To match this <p dependence, we write the matter source term as a similar Fourier decomposition 

T = T m (r,t,6) e imfl >. 
With a modicum of algebra one can rewrite the Tcukolsky equation as, 
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This is possible with the following definitions 
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where S 2 = (r 2 + a 2 ) 2 - Aa 2 sin 2 9. 

One now introduces an auxiliary field II that converts the Teukolsky equation to a set of 2 coupled PDEs which 
are first order in space and time 

~dt + dr* = ( ~ dr*^dr* + shi9~d8 [«* 6 B0 ) ~ ~ ™ ~ ^ ^ 

The code we are using produces stable evolutions of these equations using a two-step Lax-Wendroff method. 

As in other applications of this Teukolsky code, we impose boundary conditions at the edges of the polar computa- 
tional domain, i.e. at the black hole horizon, at the rotation axis and at the far end of the radial grid. The condition 
near the horizon is <f> = n = 0, due to the known asymptotic behavior of the fields 0,0]. At the outer boundary we 
impose usual outgoing boundary conditions. Since the wave equation has a potential, errors due to reflection from this 
boundary bounce back into the computational domain. We deal with this by making the computational domain in 
the radial direction large enough that any numerical reflections will not make it back to the point where we compute 
the waveforms and energy flux in the time allotted for the simulation. At the axis of rotation of the black hole one 
imposes either the condition $ = or dg$> = depending on the parity of the field specified by the azimuthal integer 
m. 

The initial data we take for the fields is zero. Since we are considering the Teukolsky equation with a source we 
are therefore violating the constraints of general relativity on the initial data. This would correspond to the particle 
"appearing suddenly" and therefore generates an artificial burst of radiation. To get more realistic initial data slices, 
one could try constructing initial data sets that correctly satisfy the Hamiltonian and momentum constraints with 
prescriptions like that of Bowen and York in the limit in which one black hole is very small, for instance, or 
modifications of particle limit data sets such as those of Lousto and Price [l^. These prescriptions specify the 
(gij,Kij) sets for use in standard ADM evolutions. To convert these quan tities to the set ($,c) t <I>) needed for the 
Teukolsky evolution, one can use the formulae developed in reference 20] . We plan to explore the extent of the 
benefits that this more accurate intial data specification may bring to enhance the orbital simulations in the particle 
limit case future work. 

The next step is to construct the source term for a general geodesic trajectory of the particle moving in the 
background Kerr geometry. To include the effect of radiation reaction one treats it as a slight deviation of the 
background geodesic orbital motion. 

In the time domain Teukolsky equation the main source term is given by 

T = 2p- 4 T 4 , (13) 

where 

T 4 = (A + 3 7 - 7 * + 4/i + fi*)[(5* - 2t* + 2a)T nm , 
- (A + 2 7 - 2 7 * + fi*)T m * m *} + (S* -r*+l3*+3a + 4tt) 

x [(A + 2 7 + 2 M *)T„ m . - (S* -t* + 2/3* + 2a)T nn ]. (14) 
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All the differential operators and Newman-Penrose coefficients, as well as the choice of Kinnersley tetrad vectors 
follow the notation of \L7\ . By choosing their values in Boyer-Lindquist coordinates and expanding everything we can 
end up with an explicit but quite complicated form of the source term that is amenable for numerical calculations |2ll| . 

We obtain the r(t), 0(t), and (j>(t) by discretizing and integrating the geodesic equations in the Kerr background, 
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The source term is singular in those points where the particle is located at each timestep. We treat this |2g] by 
using the approximation in which we substitute the delta function by a very narrow Gaussian function (whose width 
is only of a few gridpoints) : 
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(19) 



One must show that this substitution is adequate by showing convergence of the obtained waveform as a — > 0, and 
one must take special care in normalizing the mass correctly as the particle moves in the background Kerr spacetimc 
so that integration over the whole Gaussian density profile gives a constant mass throughout the whole trajectory 
and to satisfy V^T^ = 0. This translates in that the Gaussian must approximately integrates to one on the spatial 
manifold. To achieve this with the correct three-metric of the background black hole spacetime, we normalize the 
mass-density of the particle by the factor, 



N = 



(20) 



where 7^) is the 3-metric of flat space and is the 3-metric of a slice of constant Boyer-Lindquist time of the 
Kerr spacetime. We then normalize the particle mass as it moves thru the space by multiplying the quantity /1 in the 
source term expression by this factor N. We handle the Dirac deltas in the radial and angular direction through the 
Gaussian approximation, whereas the 5{<f> — 4>(t)) function we handled analytically. The normalizations presented are 
common in the SPH literature 1231. 



III. NUMERICAL RESULTS 



We now want to calculate the energy flux emitted in the form of gravitational waves from the binary that an 
orbiting observatory such as LISA would be able to detect. In this paper we focus on calculating fluxes from circular 
equatorial orbits without the effects of radiation reaction mainly to test the viability of this approach and to show 
stability, convergence, and agreement with previous results of high accuracy frequency domain calculations!^. In 
the next section we will explain the advantages of working in the time domain for future calculations of eccentric and 
inclined particle orbits. 

The code monitors the value of tp at a set point in the grid far from the horizon and computes the energy and 
angular momentum fluxes in gravitational radiation according to the formulas |24l l25j : 
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FIG. 1: Analysis of convergence factor C calculated from evolved data at a point at a radius of 100 M in the equatorial plane 
of the black hole. Results show that the code has second order convergence (factor=4). Oscillations in the convergence factor 
for the fields measured at a point is common in systems where the solution itself oscillates. One can get a smoother plot by 
showing convergence of quantities in L2 norm, for example. 



The code produces output files listing (a) the real and imaginary parts of the ip field at representative locations of 
the 2-D grid, and (b) the computed energy and angular momentum fluxes for different values of m. And all this is 
done periodically at preselected times along the length of the simulation for a circular orbital geodesic and from here 
the fluxes are postprocessed. 

We made standard convergence tests by assuming that the numerical values included errors that decrease as a 
power of the grid spacing h. Specifically we assume that the value of the Teukolsky function $ at any point in the 
computational domain goes like 

- ^exact + Ch a + Q(h a+1 ). (23) 

Grids typically used in this work are — 40M < r* < 460Af and < 6j < ir and i ~ 6400 and j ~ 60. Therefore, if we 
run at various resolutions (L: 2000X20, M: 4000X40, H: 8000X80) we can get a quantity C that measures the order of 
convergence a of the code. This is defined as, 

<£>r — 

C = 7T-> ( 24 ) 

and results for it are depicted in figure 1 where we see second order convergence of the evolution code. We also verified 
that the results are independent of the Gaussian width a for both the radial and angular Gaussian distributions used 
to approximate the location of the particle in the discrete spatial grid (see figure 5). 

Since equatorial orbits have been treated with great precision in the frequency domain method, it is important 
to show that our method gives comparable results if one looks at quantities like the average gravitational energy 
flux at large distances from the horizon. Finn and Thorne |23j | published a comprehensive work that tabulates and 
includes all previous results for high-precision gravitational wave flux and wave amplitudes coming from many kinds of 
representative compact objects in equatorial circular orbits around massive central black holes. In Table ITT1 we present 
representative comparisons with the results of formulas (3.8) and (3.10) of that paper using the precise relativistic 
corrections included in Tables III - VI that they also provide. 



IV. DISCUSSION 



We have shown that it is possible to integrate the Teukolsky equation in the time domain for the problem of a particle 
orbiting around a black hole with the particle was represented as a Gaussian and obtain results reasonably close to 
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in mode 


Tihig doni&in Gncrgy flux 


Finn- Thorn© flux 


i 


2.5 x 1(T 1() 


2.2 x 10" 10 


2 


1.0 x 10~ 07 


i A in — 07 

1.4 x 10 


3 


3.9 x 10" 08 


3.1 x 10" 08 


4 


1.0 x 10~ 08 


8.0 x 10" 09 



TABLE II: Comparisons of gravitational wave energy fluxes detected at infinity using the frequency domain solution of the 
Teukolsky-Sasaki-Nakamura equation as calculated in Ref. 23] with the results measured numerically at r = 100A/ while 

evolving the Teukolsky equation in the time domain under the "particle-as -a-gaussian-distribution" approximation discussed 
here. The particle here is in a circular, equatorial orbit of radius 2.0ri SCO about a black hole with Kerr parameter a/M — 0.9 
and fi/M = 0.01. 




time 



FIG. 2: Representative data showing the evolution of the real part of the Teukolsky function $ for the m = 3 mode of a particle 
orbiting a rapidly rotating Kerr hole of a = 0.9M at a radius of 2ri SCO {ri SCO = 2.32M) at an observation point r/M = 100, 
= 7r/2. 



those of the frequency-domain calculations. The time domain approach is expected to become more useful for the 
case of more complicated orbits. We plan in future papers to discuss how to how to use better approximations to the 
Gaussian like those of reference [2^|, and how to incorporate radiation-reaction effects. Also, the small discrepancies 
with the Thorne-Finn results need to be further investigated to define well a "window of reliability" for the results 
of the code. Another element to be incorporated is to take into account the radiation flux falling into the black hole, 
since this will also influence trajectories corrected by radiation-reaction. The first approach will be to correct the 
orbit trajectory using energy-balance arguments, as for instance in [b4ll27|. Further possibilities include incorporating 
proposed radiation reaction forces directly into the integration of the orbit. 
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FIG. 4: Teukolsky function insensitivity to the radial and angular width of the gaussian distribution that models the particle 
in orbit. The parameters for this evolution were the same as used for the other figures. The solid curve corresponds to radial 
and angular widths of 0.75M and 0.47 rad respectively, while dotted curve corresponds to widths of 0.25M and 0.16 rad. 
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traditional sense when one makes the Gaussian narrower, since at some point it will become too narrow for the finite resolution 
of the code to model it. What can be expected, and is shown in the figure is that the results tend to a constant value for a 
while when the Gaussian is made narrower, and shortly before it becomes too narrow for the code to resolve. 
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An alternate way to treat the polar angle dependence of the source term is to expand it using spin-weighted spherical 
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